| Descargar notebook.

5. Compilación de Resultados y Validación

5.1. Cargar las bibliotecas requeridas

El siguiente código carga las bibliotecas requeridas.

import numpy as np                                            # ndarrays for gridded data
import pandas as pd                                           # DataFrames for tabular data

import geostatspy.GSLIB as GSLIB                              # GSLIB utilities, visualization and wrapper
import geostatspy.geostats as geostats                        # GSLIB methods convert to Python
import geostatspy
print('GeostatsPy version: ' + str(geostatspy.__version__))

import os                                                     # set working directory, run executables

from tqdm import tqdm                                         # suppress the status bar
from functools import partialmethod
tqdm.__init__ = partialmethod(tqdm.__init__, disable=True)

ignore_warnings = True                                        # ignore warnings?

import matplotlib.pyplot as plt                               # for plotting
import matplotlib as mpl                                      # custom colorbar
import matplotlib.ticker as mticker                           # custom colorpar ticks
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
plt.rc('axes', axisbelow=True)                                # plot all grids below the plot elements
if ignore_warnings == True:
    import warnings
    warnings.filterwarnings('ignore')

import os
import matplotlib.image as mpimg
import subprocess
import time
from matplotlib.offsetbox import AnchoredText
import matplotlib.gridspec as gridspec

import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

# trimmed inverted rainbow colormap hsv removing magenta tones
cmap = mpl.colors.ListedColormap(mpl.cm.hsv(np.linspace(0.0, 0.70, 256))).reversed()

GeostatsPy version: 0.0.78

def GSLIB2ndarray_3D(data_file, kcol,nreal, nx, ny, nz):
    if nz > 1 and ny > 1:
        array = np.ndarray(shape = (nreal, nz, ny, nx), dtype=object, order="F")
    elif ny > 1:
        array = np.ndarray(shape=(nreal, ny, nx), dtype=object, order="F")
    else:
        array = np.ndarray(nreal, nx, dtype=object)

    with open(data_file) as f:
        head = [next(f) for _ in range(2)]  # read first two lines
        line2 = head[1].split()
        ncol = int(line2[0])  # get the number of columns

        for icol in range(ncol):  # read over the column names
            head = next(f)
            if icol == kcol:
                col_name = head.split()[0]
        for ineal in range(nreal):
            if nz > 1 and ny > 1:
                for iz in range(nz):
                    for iy in range(ny):
                        for ix in range(nx):
                            head = next(f)
                            array[ineal][nz - 1 - iz][ny - 1 - iy][ix] = head.split()[kcol]
            elif ny > 1:
                for iy in range(ny):
                    for ix in range(0, nx):
                        head = next(f)
                        array[ineal][ny - 1 - iy][ix] = head.split()[kcol]
            else:
                for ix in range(nx):
                    head = next(f)
                    array[ineal][ix] = head.split()[kcol]
    return array, col_name

Nuestro punto de partida sera la base de datos compositada a 6 metros y simplificada pero completa, es decir, sin filtrar aquellas muestras por UG. Esta base de datos se utilizara para validar la estimación, verificando que las estadísticas de los bloques sean coherentes con las estadísticas de las muestras.

## Cargar data simplificada de la sección anterior

import pandas as pd
import numpy as np

# Cargar base de datos en pandas y aproximar a dos decimales
df = pd.read_csv('data.csv', sep=',', encoding='latin1').round(2)
df['UG'] = df['UG'].astype(int)  # Asegurarse de que la columna 'UG' sea de tipo entero

df # mostrar la data

	X	Y	Z	cu_pct	UG
0	472187.34	6925805.69	4211.77	0.02	2
1	472187.74	6925806.53	4205.84	0.01	2
2	472188.09	6925807.38	4199.92	0.01	2
3	472188.47	6925808.28	4194.00	0.01	2
4	472188.88	6925809.22	4188.08	0.01	2
...	...	...	...	...	...
16877	471922.69	6925391.61	4022.52	0.09	4
16878	471924.63	6925394.53	4017.65	0.11	4
16879	471926.59	6925397.45	4012.78	0.18	4
16880	471928.55	6925400.36	4007.92	0.15	4
16881	471930.53	6925403.26	4003.05	0.21	4

16882 rows × 5 columns

Procedemos a compilar algunas estadísticas básicas de las muestras por UG.

# Estadísticas por UG (Number, Max, Min, Mean, Variance)
# Asegúrese de haber ejecutado previamente la celda que carga 'df' (data.csv)
import pandas as pd

# lista de UGs a reportar (ajustar si hay más o diferentes!!!)
ug_list = [2, 3, 4]

stats = {}
for ug in ug_list:
    serie = df.loc[df['UG'] == ug, 'cu_pct']
    stats[f'UG {ug}'] = serie.agg(['count', 'max', 'min', 'mean', 'var'])

# estadisticas globales (todas las muestras)
global_stats = df['cu_pct'].agg(['count', 'max', 'min', 'mean', 'var'])

# montar DataFrame final con el orden solicitado - AJUSTAR SI HAY MÁS UGS!!!
df_stats = pd.DataFrame({
    'UG 2': stats['UG 2'],
    'UG 3': stats['UG 3'],
    'UG 4': stats['UG 4'],
    'Global': global_stats
})

# Renombrar índices a los labels solicitados
df_stats.index = ['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']

# Mostrar con 2 decimales (ajustar si lo desea)
from IPython.display import display
display(df_stats.round(2))

# Opcional: exportar a csv para revisión fuera del notebook y ajustar decimales
# df_stats.round(6).to_csv('./03_resultados/summary_stats_UGs.csv')

	UG 2	UG 3	UG 4	Global
Número	6203.00	2048.00	8631.00	16882.00
Máximo	2.03	1.64	6.72	6.72
Mínimo	0.00	0.00	0.00	0.00
Media	0.06	0.13	0.25	0.17
Varianza	0.01	0.02	0.05	0.04

5.2. Cargar modelos de UGs y Leyes

El siguiente código carga nuestras estimaciones previas, tanto de UGs como de leyes.

xmin = 471116.24; xmax = 473023.89                  # rango de valores x
ymin = 6924721.57; ymax = 6926164.63                        # rango de valores y
zmin = 2950.02; zmax = 4445.78                              # rango de valores z

xsiz = 15; ysiz = 15; zsiz = 15                             # tamaño de celda
nx = 128; ny = 97; nz = 100                         # número de celdas

tmin = -999; tmax = 999;                                      # límites de recorte de datos

# calculo automático del centro de las capas y cortes para plotear, no modificar

icx = int(nx/2); cx = xmin + icx*xsiz
icy = int(ny/2); cy = ymin + icy*ysiz
icz = int(nz/2); cz = zmin + icz*zsiz

# Cargar modelos de UGs
modelo_ug = GSLIB2ndarray_3D("./03_resultados/modelo_ug.dat ", 0, 1, nx, ny, nz)[0][0].astype(int)

# Cargar modelos de leyes

estCU_UG2 = GSLIB2ndarray_3D("./03_resultados/modelo_leyes_UG2.dat ", 0, 1, nx, ny, nz)[0][0].astype(float)
categoria_UG2 = GSLIB2ndarray_3D("./03_resultados/modelo_categorias_UG2.dat ", 0, 1, nx, ny, nz)[0][0]

estCU_UG3 = GSLIB2ndarray_3D("./03_resultados/modelo_leyes_UG3.dat ", 0, 1, nx, ny, nz)[0][0].astype(float)
categoria_UG3 = GSLIB2ndarray_3D("./03_resultados/modelo_categorias_UG3.dat ", 0, 1, nx, ny, nz)[0][0]

estCU_UG4 = GSLIB2ndarray_3D("./03_resultados/modelo_leyes_UG4.dat ", 0, 1, nx, ny, nz)[0][0].astype(float)
categoria_UG4 = GSLIB2ndarray_3D("./03_resultados/modelo_categorias_UG4.dat ", 0, 1, nx, ny, nz)[0][0]

## Combinar modelos de leyes en un solo modelo final de acuerdo a las UGs
modelo_leyes_final = np.full((nz, ny, nx), np.nan)  # inicializar con NaN
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            ug_value = modelo_ug[iz, iy, ix]
            if ug_value == 2:
                modelo_leyes_final[iz, iy, ix] = estCU_UG2[iz, iy, ix]
            elif ug_value == 3:
                modelo_leyes_final[iz, iy, ix] = estCU_UG3[iz, iy, ix]
            elif ug_value == 4:
                modelo_leyes_final[iz, iy, ix] = estCU_UG4[iz, iy, ix]

## Combinar modelos de categorías en un solo modelo final de acuerdo a las UGs
modelo_categorias_final = np.full((nz, ny, nx), "Sin_Categoria", dtype=object)  # inicializar con "Sin_Categoria"
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            ug_value = modelo_ug[iz, iy, ix]
            if ug_value == 2:
                modelo_categorias_final[iz, iy, ix] = categoria_UG2[iz, iy, ix]
            elif ug_value == 3:
                modelo_categorias_final[iz, iy, ix] = categoria_UG3[iz, iy, ix]
            elif ug_value == 4:
                modelo_categorias_final[iz, iy, ix] = categoria_UG4[iz, iy, ix]

continuous = 'cu_pct'

x_coords = np.linspace(xmin + xsiz/2, xmax - xsiz/2, nx)
y_coords = np.linspace(ymin + ysiz/2, ymax - ysiz/2, ny)
z_coords = np.linspace(zmin + zsiz/2, zmax - zsiz/2, nz)

import plotly.io as pio
# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

fig = go.Figure(data=[go.Scatter3d(
    x=df['X'],
    y=df['Y'],
    z=df['Z'],
    marker=dict(color=df[continuous],
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=df[continuous].quantile(0.95),
                size=2.0, # Esta opcion controla el tamaño de los datos
                colorbar=dict(
                    title='Cu [%]',
                    thickness=20,
                    # Add a border to the colorbar
                    outlinecolor='black',  # Color of the border
                    outlinewidth=2,       # Width of the border
                    bordercolor='white',   # Background border color (if applicable)
                    borderwidth=1         # Background border width
                )
                ),
    mode='markers',
    opacity=1
)])

# otras opciones para controlar el color de los elementos de la figura, borrar hovermode para desactivar etiquetas al pasar el mouse
fig.update_layout(    title_text='Estimación de Cu - Global',
scene = dict(
                    xaxis = dict(
                        title='Este',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    yaxis = dict(
                        title='Norte',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white"),
                    zaxis = dict(
                        title='Elevación',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    ),
                    font_family="Times New Roman",
                    font_color="black", hovermode=False
                  )

# ocultar la leyenda/autonombre "trace 0" junto al colorbar
fig.update_layout(showlegend=False)

# definir cortes en x, y, z

icz = 60  # índice de corte en z (puede ajustarse), va entre 0 y nz-1

# Agregar cortes del modelo de UGs con la misma escala de colores de plotly
x_slice = np.flip(np.flipud(modelo_leyes_final[:, :, icx]), axis=1)
y_slice = np.flip(modelo_leyes_final[:, icy, :], axis=0)
z_slice = np.flip(modelo_leyes_final[nz-icz, :, :], axis=0)

# X-slice (plano constante x)
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])

surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Y-slice (plano constante y)
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])

surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None,
)

# Z-slice (plano constante z)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])

surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Añadir superficies al figure (manteniendo la leyenda y colores consistentes)
fig.add_trace(surf_x)
fig.add_trace(surf_y)
fig.add_trace(surf_z)

fig.show()

import plotly.graph_objects as go
import plotly.express as px
from plotly.subplots import make_subplots

ROWS = 3
COLS = 2
# Crear listado de bancos
step_size = 15
max_benches = ROWS*COLS
lims = np.arange(4250,4250-step_size*max_benches,-step_size)

continuous = 'cu_pct'

# Initialize figure with subplots
# Initialize figure with subplots
fig = make_subplots(
    rows=ROWS, cols=COLS, horizontal_spacing=0.2, vertical_spacing=0.1,
    subplot_titles=[f"Cu [%], banco {lims[i]+step_size/2}" for i in range(len(lims))]
)

for i, lim in enumerate(lims):
    data_temp = df.loc[(df['Z'] > lim) & (df['Z'] <= lim + step_size)]
    # Calculate row and column indices for the subplot (1-based indexing)
    row = i // COLS + 1
    col = i % COLS + 1

    # agregar Z-slice del modelo_leyes_final correspondiente al banco (plano constante z)
    z_mid = lim + step_size/2.0
    iz = int(np.abs(z_coords - z_mid).argmin())  # índice de la capa más cercana al centro del banco

    # extraer y orientar la rebanada (debe quedar shape (ny, nx) con y_coords como filas)
    z_slice2d = np.flip(modelo_leyes_final[nz-iz, :, :].astype(float), axis=0)

    # enmascarar valores inválidos (-999) para que sean transparentes en el heatmap
    z_slice2d[z_slice2d <= -998] = np.nan

    heat = go.Heatmap(
        x=x_coords,
        y=y_coords,
        z=z_slice2d,
        colorscale=px.colors.sequential.Rainbow[1:],
        zmin=0.0,
        zmax=df[continuous].quantile(0.95),
        showscale=False,
        opacity=0.7,
        hoverinfo='skip',
        name=f'Z-slice {z_coords[iz]:.1f}'
    )

    # añadir primero el heatmap (fondo) y luego el scatter (muestras) conserva la visibilidad
    fig.add_trace(heat, row=row, col=col)

    fig.add_trace(
        go.Scatter(
            x=data_temp['X'],
            y=data_temp['Y'],
            marker=dict(
                color=data_temp[continuous],
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=df[continuous].quantile(0.95),
                size=5.0,
                line=dict(
                    color='black',
                    width=0.5
                ),
            ),
            mode='markers',
            opacity=1,
            showlegend=False
        ),
        row=row, col=col
    )



    # Update xaxis properties
    fig.update_xaxes(title_text="Este", range=[df['X'].quantile(0.015), df['X'].quantile(0.97)], gridcolor='rgba(0,0,0,0.2)', row=row, col=col)

    # Update yaxis properties
    fig.update_yaxes(title_text="Norte", range=[df['Y'].quantile(0.03), df['Y'].quantile(0.99)], gridcolor='rgba(0,0,0,0.2)', row=row, col=col)

fig.update_annotations(font_size=12)
# otras opciones para controlar el color de los elementos de la figura
fig.update_layout(font_family="Times New Roman",font_size=8,
                    font_color="black", plot_bgcolor='white',   autosize = False,
    width = 600,
    height = 900, hovermode=False)


fig.show()

from plotly.subplots import make_subplots

## truncate Z variable in the df axis every x meters

df_swath = df.copy()

dz = 50 # truncar cada 50 metros
df_swath['Z'] = df_swath['Z'].where(df_swath['Z'] % dz == 0, other=df_swath['Z'] - df_swath['Z'] % dz)

## Plot truncated z values against filtered average cut > 0 percentage with a line plot and marginal counter histogram
import plotly.express as px
import plotly.graph_objects as go

# prepare aggregated line data
agg = df_swath[df_swath['cu_pct'] >= 0].groupby('Z').agg(
    mean_cut=('cu_pct', 'mean'),
    count=('cu_pct', 'count')
).reset_index()

# create subplot: top = line, bottom = marginal histogram (shared x)
# enable a secondary y-axis for the bottom subplot so we can plot the bar on the right
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
                    vertical_spacing=0.05,
                    row_heights=[0.7, 0.3],
                    specs=[[{"type": "xy"}],
                           [{"type": "xy", "secondary_y": True}]])

# line trace (top)
fig.add_trace(
    go.Scatter(
        x=agg['Z'],
        y=agg['mean_cut'],
        mode='lines+markers',
        name='Ley Promedio Cu [%] en muestras',
        line=dict(color='blue'),
        marker=dict(size=6, color='blue')
    ),
    row=1, col=1
)

## add z vs filtered line
# line trace (top)
fig.add_trace(
    go.Scatter(
        x=z_coords,
        y=[np.mean(np.mean(modelo_leyes_final[nz-iz-1, :, :][modelo_leyes_final[nz-iz-1, :, :] > 0])) for iz in range(nz)],
        mode='lines',
        name='Ley Promedio Cu (%) en modelo',
        line=dict(color='red')
    ),
    row=1, col=1
)

# marginal histogram (bottom)
fig.add_trace(
    go.Histogram(
        x=df_swath[df_swath['cu_pct'] >= 0]['Z'],
        nbinsx=30,
        marker=dict(color='blue'),
        opacity=0.3,
        name='Conteo muestras'
    ),
    row=2, col=1
)

# secondary axis: count of valid cells (>0) per model z-layer
valid_counts = [
    int(np.sum((modelo_leyes_final[nz - iz - 1, :, :] > 0) & np.isfinite(modelo_leyes_final[nz - iz - 1, :, :])))
    for iz in range(nz)
]

bar = go.Bar(
    x=z_coords,
    y=valid_counts,
    name='N° celdas válidas (modelo)',
    marker=dict(color='red'),
    opacity=0.3
)
## add bar y-axis on the right, row=2, col=1.

# add the bar to the bottom subplot using the secondary y-axis
fig.add_trace(bar, row=2, col=1, secondary_y=True)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
    title_text='Celdas válidas (modelo)',
    secondary_y=True,
    row=2, col=1
)

# layout and axis labels
fig.update_xaxes(title_text='Elevación (m)', row=2, col=1)
fig.update_yaxes(title_text='Ley Promedio Cu (%)', row=1, col=1)
fig.update_yaxes(title_text='Conteo muestras', row=2, col=1)
fig.update_yaxes(title_text='Celdas válidas (modelo)', row=2, col=1, secondary_y=True)

fig.update_layout(
    title='Ley Promedio de Cu vs Elevación (m)',
    font_family="Times New Roman",
    font_size=12,
    font_color="black",
    plot_bgcolor='white',
    hovermode=False,
    showlegend=True,
    bargap=0.1,
    height=600,
    width=900
)
fig.show()

## esta vez para UG 4

from plotly.subplots import make_subplots

## truncate Z variable in the df axis every x meters

df_swath = df[df['UG'] == 4].copy()

dz = 50 # truncar cada 50 metros
df_swath['Z'] = df_swath['Z'].where(df_swath['Z'] % dz == 0, other=df_swath['Z'] - df_swath['Z'] % dz)

## Plot truncated z values against filtered average cut > 0 percentage with a line plot and marginal counter histogram
import plotly.express as px
import plotly.graph_objects as go

# prepare aggregated line data
agg = df_swath[df_swath['cu_pct'] >= 0].groupby('Z').agg(
    mean_cut=('cu_pct', 'mean'),
    count=('cu_pct', 'count')
).reset_index()

# create subplot: top = line, bottom = marginal histogram (shared x)
# enable a secondary y-axis for the bottom subplot so we can plot the bar on the right
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
                    vertical_spacing=0.05,
                    row_heights=[0.7, 0.3],
                    specs=[[{"type": "xy"}],
                           [{"type": "xy", "secondary_y": True}]])

# line trace (top)
fig.add_trace(
    go.Scatter(
        x=agg['Z'],
        y=agg['mean_cut'],
        mode='lines+markers',
        name='Ley Promedio Cu [%] en muestras',
        line=dict(color='blue'),
        marker=dict(size=6, color='blue')
    ),
    row=1, col=1
)

modelo_leyes_final_UG4 = np.full((nz, ny, nx), np.nan)  # inicializar con NaN
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            if modelo_ug[iz, iy, ix] == 4:
                modelo_leyes_final_UG4[iz, iy, ix] = estCU_UG4[iz, iy, ix]


## add z vs filtered line
# line trace (top)
fig.add_trace(
    go.Scatter(
        x=z_coords,
        y=[np.mean(np.mean(modelo_leyes_final_UG4[nz-iz-1, :, :][modelo_leyes_final_UG4[nz-iz-1, :, :] > 0])) for iz in range(nz)],
        mode='lines',
        name='Ley Promedio Cu (%) en modelo',
        line=dict(color='red')
    ),
    row=1, col=1
)

# marginal histogram (bottom)
fig.add_trace(
    go.Histogram(
        x=df_swath[df_swath['cu_pct'] >= 0]['Z'],
        nbinsx=30,
        marker=dict(color='blue'),
        opacity=0.3,
        name='Conteo muestras'
    ),
    row=2, col=1
)

# secondary axis: count of valid cells (>0) per model z-layer
valid_counts = [
    int(np.sum((modelo_leyes_final_UG4[nz - iz - 1, :, :] > 0) & np.isfinite(modelo_leyes_final_UG4[nz - iz - 1, :, :])))
    for iz in range(nz)
]

bar = go.Bar(
    x=z_coords,
    y=valid_counts,
    name='N° celdas válidas (modelo)',
    marker=dict(color='red'),
    opacity=0.3
)
## add bar y-axis on the right, row=2, col=1.

# add the bar to the bottom subplot using the secondary y-axis
fig.add_trace(bar, row=2, col=1, secondary_y=True)
# Update 2nd y-axis properties for the bar chart
fig.update_yaxes(
    title_text='Celdas válidas (modelo)',
    secondary_y=True,
    row=2, col=1
)

# layout and axis labels
fig.update_xaxes(title_text='Elevación (m)', row=2, col=1)
fig.update_yaxes(title_text='Ley Promedio Cu (%)', row=1, col=1)
fig.update_yaxes(title_text='Conteo muestras', row=2, col=1)
fig.update_yaxes(title_text='Celdas válidas (modelo)', row=2, col=1, secondary_y=True)

fig.update_layout(
    title='Ley Promedio de Cu vs Elevación (m) - UG 4',
    font_family="Times New Roman",
    font_size=12,
    font_color="black",
    plot_bgcolor='white',
    hovermode=False,
    showlegend=True,
    bargap=0.1,
    height=600,
    width=900
)
fig.show()

from IPython.display import display

# entregar estadísticas de modelo_leyes_final según UG (Number, Max, Min, Mean, Variance)

ugs = [2, 3, 4]

results = {}
for ug in ugs:
    mask = (modelo_ug == ug) & (modelo_leyes_final > 0)
    vals = modelo_leyes_final[mask]
    cnt = vals.size
    if cnt > 0:
        vals = vals.astype(float)
        results[f'UG {ug}'] = {
            'Número': int(cnt),
            'Máximo': float(np.nanmax(vals)),
            'Mínimo': float(np.nanmin(vals)),
            'Media': float(np.nanmean(vals)),
            'Varianza': float(np.nanvar(vals))
        }
    else:
        results[f'UG {ug}'] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}

# Estadísticas globales (todas las celdas válidas)
mask_all = (modelo_leyes_final > 0) & np.isfinite(modelo_leyes_final)
vals_all = modelo_leyes_final[mask_all]
cnt_all = vals_all.size
if cnt_all > 0:
    results['Global'] = {
        'Número': int(cnt_all),
        'Máximo': float(np.nanmax(vals_all)),
        'Mínimo': float(np.nanmin(vals_all)),
        'Media': float(np.nanmean(vals_all)),
        'Varianza': float(np.nanvar(vals_all))
    }
else:
    results['Global'] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}

# Montar DataFrame y mostrar
df_stats_model = pd.DataFrame(results).T[['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']]
display(df_stats_model.round(2))

	Número	Máximo	Mínimo	Media	Varianza
UG 2	125529.0	0.73	0.00	0.05	0.00
UG 3	20092.0	0.70	0.01	0.12	0.01
UG 4	177774.0	2.46	0.00	0.22	0.01
Global	323395.0	2.46	0.00	0.15	0.02

# Ensure renderer is appropriate for your environment (notebook / vscode)
pio.renderers.default = 'notebook'

# --- Input assumptions (these must exist in the environment) ---
# df: pandas DataFrame with columns 'X','Y','Z'
# categoria_UG4: numpy array with shape (nz, ny, nx) containing strings: 'Inferido','Indicado','Medido'
# x_coords, y_coords, z_coords: 1D coordinate arrays for grid cell centers (lengths nx, ny, nz)
# icx, icy, icz: integer indices for the central slices to display
# nz, ny, nx: grid dimensions

# Map categories to numeric codes
cat_to_num = {'Inferido': 0, 'Indicado': 1, 'Medido': 2}
cat_num = np.vectorize(lambda x: cat_to_num.get(x, np.nan))(modelo_categorias_final)
cat_num = cat_num.astype(float)  # allow NaNs

# Create slices (same convention used in the notebook)
x_slice = np.flip(np.flipud(cat_num[:, :, icx]), axis=1)   # X-slice (nz, ny)
y_slice = np.flip(cat_num[:, icy, :], axis=0)             # Y-slice (nz, nx)
z_slice = np.flip(cat_num[nz-icz, :, :], axis=0)          # Z-slice (ny, nx)

# Discrete colors: Inferido (red), Indicado (yellow), Medido (green)
# We use normalized stops [0.0, 0.5, 1.0] so numeric codes 0,1,2 map exactly to the three colors
colorscale = [
    [0.0, '#d73027'],   # Inferido -> 0
    [0.5, "#fff700"],   # Indicado -> 1 (normalized to 0.5)
    [1.0, '#1a9850']    # Medido -> 2 (normalized to 1.0)
]

# Build surface for X-slice
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
## nan si no hay datos en la celda
XX[np.isnan(x_slice)] = np.nan
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Build surface for Y-slice
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
## nan si no hay datos en la celda
YY[np.isnan(y_slice)] = np.nan
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Build surface for Z-slice (we enable a colorbar here with categorical ticks)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
## nan si no hay datos en la celda
ZZ[np.isnan(z_slice)] = np.nan
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=True,
    colorbar=dict(
        title='Categoría',
        tickvals=[0, 1, 2],
        ticktext=['Inferido', 'Indicado', 'Medido'],
        thickness=20
    ),
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Add sample points as neutral black markers for reference
scatter_samples = go.Scatter3d(
    x=df['X'], y=df['Y'], z=df['Z'],
    mode='markers',
    marker=dict(color='black', size=2, opacity=0.8),
    name='Muestras'
)

# Compose figure
fig = go.Figure(data=[scatter_samples, surf_x, surf_y, surf_z])
fig.update_layout(
    title='Categoría (Inferido=rojo, Indicado=amarillo, Medido=verde)',
    scene=dict(
        xaxis=dict(title='Este'),
        yaxis=dict(title='Norte'),
        zaxis=dict(title='Elevación')
    ),
    width=900,
    height=720,
    font=dict(family='Times New Roman', size=12)
)

fig.show()

# Entregar para cada UG y cada categoría las estadísticas (Número, Max, Min, Mean, Variance) de modelo_leyes_final
from IPython.display import display

ugs = [2, 3, 4]
categories = ['Inferido', 'Indicado', 'Medido']
results = {}
for ug in ugs:
    for cat in categories:
        mask = (modelo_ug == ug) & (modelo_categorias_final == cat) & (modelo_leyes_final > 0)
        vals = modelo_leyes_final[mask]
        cnt = vals.size
        key = f'UG {ug} - {cat}'
        if cnt > 0:
            vals = vals.astype(float)
            results[key] = {
                'Número': int(cnt),
                'Máximo': float(np.nanmax(vals)),
                'Mínimo': float(np.nanmin(vals)),
                'Media': float(np.nanmean(vals)),
                'Varianza': float(np.nanvar(vals))
            }
        else:
            results[key] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}
# Montar DataFrame y mostrar
df_stats_cat = pd.DataFrame(results).T[['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']]
display(df_stats_cat.round(2))

	Número	Máximo	Mínimo	Media	Varianza
UG 2 - Inferido	64765.0	0.73	0.00	0.05	0.00
UG 2 - Indicado	37574.0	0.72	0.00	0.06	0.00
UG 2 - Medido	23190.0	0.71	0.00	0.06	0.01
UG 3 - Inferido	6664.0	0.43	0.01	0.10	0.00
UG 3 - Indicado	7136.0	0.46	0.01	0.12	0.01
UG 3 - Medido	6292.0	0.70	0.01	0.13	0.01
UG 4 - Inferido	88448.0	1.47	0.03	0.21	0.01
UG 4 - Indicado	69085.0	1.50	0.01	0.23	0.01
UG 4 - Medido	20241.0	2.46	0.00	0.25	0.02

# Entregar para cada UG y cada categorías global = med+ind+inf y demostrado = med + ind las estadísticas (Número, Max, Min, Mean, Variance) de modelo_leyes_final
from IPython.display import display

ugs = [2, 3, 4]
categories = [['Inferido', 'Indicado', 'Medido'], ['Indicado', 'Medido']]
category_labels = ['Global', 'Demostrado']
results = {}

for ug in ugs:
    for cat, label in zip(categories, category_labels):
        mask = (modelo_ug == ug) & np.isin(modelo_categorias_final, cat) & (modelo_leyes_final > 0) & np.isfinite(modelo_leyes_final)
        vals = modelo_leyes_final[mask]
        cnt = vals.size
        key = f'UG {ug} - {label}'
        if cnt > 0:
            vals = vals.astype(float)
            results[key] = {
                'Número': int(cnt),
                'Máximo': float(np.nanmax(vals)),
                'Mínimo': float(np.nanmin(vals)),
                'Media': float(np.nanmean(vals)),
                'Varianza': float(np.nanvar(vals))
            }
        else:
            results[key] = {'Número': 0, 'Máximo': np.nan, 'Mínimo': np.nan, 'Media': np.nan, 'Varianza': np.nan}
# Montar DataFrame y mostrar
df_stats_cat = pd.DataFrame(results).T[['Número', 'Máximo', 'Mínimo', 'Media', 'Varianza']]
display(df_stats_cat.round(2))

	Número	Máximo	Mínimo	Media	Varianza
UG 2 - Global	125529.0	0.73	0.00	0.05	0.00
UG 2 - Demostrado	60764.0	0.72	0.00	0.06	0.00
UG 3 - Global	20092.0	0.70	0.01	0.12	0.01
UG 3 - Demostrado	13428.0	0.70	0.01	0.13	0.01
UG 4 - Global	177774.0	2.46	0.00	0.22	0.01
UG 4 - Demostrado	89326.0	2.46	0.00	0.23	0.02

## Curva Tonelaje (izq.) y ley promedio (der.) in secondary axis para el modelo final de leyes overlayed in the same plot
import plotly.graph_objects as go
# Prepare data
cutoff_values = np.arange(0, df['cu_pct'].quantile(0.96), 0.01)
mtonnage = []
average_grade = []
for cutoff in cutoff_values:
    mask = (modelo_leyes_final >= cutoff) & np.isfinite(modelo_leyes_final)
    ## puede crear otra mascara mask para filtrar por categoría o UG si lo desea
    mask = np.logical_and(mask, np.logical_or(modelo_categorias_final == 'Medido', modelo_categorias_final == 'Indicado', modelo_categorias_final == 'Inferido'))  # ejemplo: solo categoría Medido
    ## mask = np.logical_and(mask, modelo_ug == 4)  # ejemplo: solo UG 4. COMBINAR AMBOS FILTROS SI ES NECESARIO

    total_tonnage = np.sum(mask) * (xsiz * ysiz * zsiz * 2.5)  # assuming 2.5 t/m³ (2.5 g/cm³)
    total_mtonnage = total_tonnage / 1e6  # convert to million tonnes
    mtonnage.append(round(total_mtonnage, 2))
    if total_tonnage > 0:
        avg_grade = np.nanmean(modelo_leyes_final[mask])
    else:
        avg_grade = 0
    average_grade.append(round(avg_grade, 2))

# Create figure with secondary y-axis
fig = go.Figure()
# Tonnage trace
fig.add_trace(
    go.Scatter(
        x=cutoff_values,
        y=mtonnage,
        name='Tonelaje',
        line=dict(color='blue'),
        yaxis='y1'
    )
)
# Average grade trace
fig.add_trace(
    go.Scatter(
        x=cutoff_values,
        y=average_grade,
        name='Ley Promedio (%)',
        line=dict(color='red'),
        yaxis='y2'
    )
)

# Update layout for dual y-axes
fig.update_layout(
    title='Curva Tonelaje y Ley Promedio vs Corte de Ley (Modelo Final de Leyes)',
    xaxis=dict(title='Corte de Ley (%)',gridcolor='rgba(0,0,0,0.2)'),
    yaxis=dict(
        title='Tonelaje (MTon)',
        tickfont=dict(color='blue',
                      ),
        side='left',
        gridcolor='rgba(0,0,0,0.2)'
    ),
    yaxis2=dict(
        title='Ley Promedio (%)',
        tickfont=dict(color='red'),
        tickformat=',.f',
        overlaying='y',
        side='right'
    ),
    legend=dict(x=0.1, y=0.9),
    font_family="Times New Roman",
    font_size=12,
    font_color="black",
    plot_bgcolor='white',
    width=800,
    height=500
)
fig.show()

6. Apéndice: Demarcación de topografía

Puede mejorar su estimación de recursos restringiendo el cálculo de ellos solo a aquellos bloques debajo de la topografía. Para esto, es necesario flagear aquellos bloques por encima/por debajo de esta, de manera de que no sean considerados.

---

6.1. Instalar librerias de manejo de CAD

!conda install -c conda-forge ezdxf scipy -y

3 channel Terms of Service accepted
Channels:
 - conda-forge
 - defaults
Platform: win-64
Collecting package metadata (repodata.json): done
Solving environment: done

# All requested packages already installed.

==> WARNING: A newer version of conda exists. <==
    current version: 25.7.0
    latest version: 25.9.1

Please update conda by running

    $ conda update -n base -c conda-forge conda

6.2. Flagear topografía en modelo de bloques vació

"""
Flag blocks below topography from a DXF surface.

Funciona con DXF que contenga 3DFACE (malla) o polylines/contours (extrae vértices).
Salida: nuevo CSV con columnas adicionales Z_topo y Below_Topo (True/False).

Requisitos: ezdxf, numpy, pandas, scipy
"""

import ezdxf
import numpy as np
import pandas as pd
from scipy.interpolate import LinearNDInterpolator
from scipy.spatial import cKDTree
import sys
import os
import geostatspy.GSLIB as GSLIB                              # GSLIB utilities, visualization and wrapper


YY, ZZ, XX = np.meshgrid(y_coords, z_coords, x_coords)             # shapes (ny, nx, nz)
## to csv
grid_points = np.column_stack([XX.ravel(), YY.ravel(), ZZ.ravel()])  # shape (N, 3) where N = nx*ny*nz
np.savetxt("grid_points.csv", grid_points, delimiter=",", header="X,Y,Z", comments="")

# ====== CONFIG ======
DXF_FILE = "topo.dxf"          # <- cambia al path de tu DXF
BLOCKS_CSV = "grid_points.csv"       # <- cambia al path de tu CSV
OUTPUT_CSV = "block_model_flagged.csv"

# Nombre de columnas del block model
X_COL = "X"
Y_COL = "Y"
Z_COL = "Z"

TOL = 1e-9   # tolerancia para comparaciones; ajústala si quieres
# ====================

def extract_surface_points_from_dxf(dxf_path):
    """
    Extrae puntos (x,y,z) desde 3DFACE, POLYLINE, LWPOLYLINE, POINT, VERTEX.
    Devuelve numpy array (N,3).
    """
    doc = ezdxf.readfile(dxf_path)
    msp = doc.modelspace()

    pts = []

    # 3DFACE y 3dface (varía según ezdxf versiones)
    try:
        faces = msp.query("3DFACE")
    except Exception:
        faces = msp.query("3dface")
    for f in faces:
        # cada face puede tener up to 4 vértices
        try:
            for v in f.vertices:
                pts.append((float(v[0]), float(v[1]), float(v[2])))
        except Exception:
            # fallback: usar atributos dxf
            for i in range(1,5):
                key = f"dxf.vtx{i}"
                if hasattr(f.dxf, f"vtx{i}"):
                    v = getattr(f.dxf, f"vtx{i}")
                    pts.append((float(v[0]), float(v[1]), float(v[2])))

    # LWPOLYLINE y POLYLINE (líneas de contorno)
    for pl in msp.query("LWPOLYLINE"):
        for p in pl.get_points():
            # get_points devuelve (x,y,bulge) o (x,y) o (x,y,z) - toleramos
            if len(p) >= 2:
                x, y = float(p[0]), float(p[1])
                z = None
                # algunos LWPOLYLINE no tienen z; si existe, ezdxf puede guardarlo en .attribs o en p[2]
                if len(p) >= 3:
                    try:
                        z = float(p[2])
                    except Exception:
                        z = None
                # si z es None, trataremos luego; por ahora, omitimos o podemos asignar 0 (mejor omitir)
                if z is not None:
                    pts.append((x, y, z))
    # POLYLINE clásicas (3D POLYLINE)
    for pl in msp.query("POLYLINE"):
        for v in pl.vertices():
            coords = v.dxf.location if hasattr(v.dxf, "location") else v
            try:
                x, y, z = float(coords[0]), float(coords[1]), float(coords[2])
                pts.append((x, y, z))
            except Exception:
                # si no hay Z explícito, ignoramos
                pass

    # POINT entities
    for p in msp.query("POINT"):
        try:
            x, y, z = p.dxf.location
            pts.append((float(x), float(y), float(z)))
        except Exception:
            pass

    pts = np.array(pts, dtype=float)
    if pts.size == 0:
        raise ValueError("No se encontraron puntos 3D en el DXF. Asegúrate de que el DXF tenga 3DFACE o polylines con Z.")
    # deduplicate (opcional)
    pts = np.unique(pts, axis=0)
    return pts

def build_interpolator(surface_pts):
    """
    Construye un interpolador LinearNDInterpolator sobre los puntos de la superficie.
    Devuelve función interp(x_array, y_array) -> z_array
    """
    xy = surface_pts[:, :2]
    z = surface_pts[:, 2]
    interp = LinearNDInterpolator(xy, z, fill_value=np.nan)
    # KDTree para nearest neighbor fill
    tree = cKDTree(xy)

    def interp_with_nn(xq, yq):
        xq = np.asarray(xq)
        yq = np.asarray(yq)
        qpts = np.column_stack([xq.ravel(), yq.ravel()])
        zq = interp(qpts)
        # interp puede devolver array con shape (N,); si NaN en puntos fuera del convex hull, rellenar con nearest
        nan_mask = np.isnan(zq)
        if np.any(nan_mask):
            dists, idx = tree.query(qpts[nan_mask], k=1)
            zq[nan_mask] = z[idx]
        return zq.reshape(xq.shape)
    return interp_with_nn

def main(dxf_file, blocks_csv, output_csv):
    print("Leyendo DXF y extrayendo puntos de superficie...")
    surface_pts = extract_surface_points_from_dxf(dxf_file)
    print(f"Puntos extraídos de la superficie: {surface_pts.shape[0]}")

    print("Construyendo interpolador de la superficie...")
    interp = build_interpolator(surface_pts)

    print("Leyendo block model CSV...")
    df = pd.read_csv(blocks_csv)
    if X_COL not in df.columns or Y_COL not in df.columns or Z_COL not in df.columns:
        print("Columnas encontradas en CSV:", df.columns.tolist())
        raise KeyError(f"El CSV debe contener columnas '{X_COL}', '{Y_COL}', '{Z_COL}'. Ajusta X_COL/Y_COL/Z_COL en el script si tus nombres son diferentes.")

    # calcular topo en cada bloque
    print("Interpolando elevación topo para cada bloque...")
    xq = df[X_COL].values
    yq = df[Y_COL].values
    z_topo = interp(xq, yq)
    df["Z_topo"] = z_topo

    # Si hay NaNs (no debería haber por nearest fill), tratarlos:
    if df["Z_topo"].isna().any():
        print("Advertencia: hay valores Z_topo NaN. Rellenando con nearest neighbor simple.")
        # fallback: nearest neighbor from surface points
        tree = cKDTree(surface_pts[:, :2])
        idx = tree.query(np.column_stack([xq, yq]), k=1)[1]
        df.loc[df["Z_topo"].isna(), "Z_topo"] = surface_pts[idx[df["Z_topo"].isna()], 2]

    # flag: block Z < topo Z  (puedes invertir el < si tus Z están ordenadas al revés)
    df["Below_Topo"] = (df[Z_COL].astype(float) < (df["Z_topo"].astype(float) - TOL)).astype(int)

    print(f"Bloques totales: {len(df)}. Bloques below topo: {df['Below_Topo'].sum()}")
    print(f"Guardando resultado en: {output_csv}")
    GSLIB.Dataframe2GSLIB(output_csv, df)
    print("Listo.")

if __name__ == "__main__":
    # si quieres pasar argumentos por línea de comando:
    # python flag_blocks_below_topo.py topography.dxf block_model.csv output.csv
    if len(sys.argv) >= 4:
        DXF_FILE = sys.argv[1]
        BLOCKS_CSV = sys.argv[2]
        OUTPUT_CSV = sys.argv[3]
    main(DXF_FILE, BLOCKS_CSV, OUTPUT_CSV)

Leyendo DXF y extrayendo puntos de superficie...
Puntos extraídos de la superficie: 14177
Construyendo interpolador de la superficie...
Leyendo block model CSV...
Interpolando elevación topo para cada bloque...
Bloques totales: 1241600. Bloques below topo: 1083384
Guardando resultado en: block_model_flagged.csv
Listo.

6.3. Leer y demarcar modelo de leyes

## leer modelo flageado
modelo_flagged = GSLIB.GSLIB2ndarray_3D("block_model_flagged.csv", 4, 1, nx, ny, nz)[0][0]

# copy modelo_leyes_final to a new array and set flagged cells to NaN
modelo_leyes_flagged = modelo_leyes_final.copy()

for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            if modelo_flagged[iz, iy, ix] == 0:  # Below_Topo == True
                modelo_leyes_flagged[iz, iy, ix] = -9999

6.4. Validación visual

continuous = 'cu_pct'

x_coords = np.linspace(xmin + xsiz/2, xmax - xsiz/2, nx)
y_coords = np.linspace(ymin + ysiz/2, ymax - ysiz/2, ny)
z_coords = np.linspace(zmin + zsiz/2, zmax - zsiz/2, nz)

import plotly.io as pio
# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

fig = go.Figure(data=[go.Scatter3d(
    x=df['X'],
    y=df['Y'],
    z=df['Z'],
    marker=dict(color=df[continuous],
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=df[continuous].quantile(0.95),
                size=2.0, # Esta opcion controla el tamaño de los datos
                colorbar=dict(
                    title='Cu [%]',
                    thickness=20,
                    # Add a border to the colorbar
                    outlinecolor='black',  # Color of the border
                    outlinewidth=2,       # Width of the border
                    bordercolor='white',   # Background border color (if applicable)
                    borderwidth=1         # Background border width
                )
                ),
    mode='markers',
    opacity=1
)])

# otras opciones para controlar el color de los elementos de la figura, borrar hovermode para desactivar etiquetas al pasar el mouse
fig.update_layout(    title_text='Estimación de Cu - Global',
scene = dict(
                    xaxis = dict(
                        title='Este',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    yaxis = dict(
                        title='Norte',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white"),
                    zaxis = dict(
                        title='Elevación',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    ),
                    font_family="Times New Roman",
                    font_color="black", hovermode=False
                  )

# ocultar la leyenda/autonombre "trace 0" junto al colorbar
fig.update_layout(showlegend=False)

# definir cortes en x, y, z

icz = 60  # índice de corte en z (puede ajustarse), va entre 0 y nz-1

# Agregar cortes del modelo de UGs con la misma escala de colores de plotly
x_slice = np.flip(np.flipud(modelo_leyes_flagged[:, :, icx]), axis=1)
y_slice = np.flip(modelo_leyes_flagged[:, icy, :], axis=0)
z_slice = np.flip(modelo_leyes_flagged[nz-icz, :, :], axis=0)

# X-slice (plano constante x)
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
# Introducir valores NaN para hacer partes de la superficie transparentes
XX[x_slice.astype(float) < -9990] = np.nan
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Y-slice (plano constante y)
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
## Introducir valores NaN para hacer partes de la superficie transparentes
YY[y_slice.astype(float) < -9990] = np.nan
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None,
)

# Z-slice (plano constante z)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
## Introducir valores NaN para hacer partes de la superficie transparentes
ZZ[z_slice.astype(float) < -9990] = np.nan
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Añadir superficies al figure (manteniendo la leyenda y colores consistentes)
fig.add_trace(surf_x)
fig.add_trace(surf_y)
fig.add_trace(surf_z)

## add 3d mesh with topography surface from dxf file for reference
topo_pts = extract_surface_points_from_dxf("topo.dxf")
topo_mesh = go.Mesh3d(
    x=topo_pts[:, 0],
    y=topo_pts[:, 1],
    z=topo_pts[:, 2],
    color='lightgrey',
    opacity=0.5,
    name='Topografía (DXF)',
    hoverinfo='skip'
)
fig.add_trace(topo_mesh)

fig.show()

## Export CSV containing id in x y z, coordinates in x y z, cu_pct, Ug , category (-99, 1 2 or 3), flagged below topo and density equal to 2.5 g/cm3 but 9999 at the all borders except at top border)
# Build a list of rows and create the DataFrame once (avoids deprecated DataFrame.append and variable shadowing)

rows = []
for iz in range(nz):
    for iy in range(ny):
        for ix in range(nx):
            x_coord = x_coords[ix]
            y_coord = y_coords[iy]
            z_coord = z_coords[iz]
            cu_pct = modelo_leyes_flagged[iz, iy, ix]
            ug = modelo_ug[iz, iy, ix]
            ## category mapping: Inferido=0, Indicado=1, Medido=2
            category_str = modelo_categorias_final[iz, iy, ix]
            if category_str == 'Inferido':
                category = 0
            elif category_str == 'Indicado':
                category = 1
            elif category_str == 'Medido':
                category = 2
            else:
                category = -99  # unknown category
            flagged_below_topo = int(modelo_flagged[iz, iy, ix] == 0)
            # density: 2.5 g/cm3 inside model, 9999 at borders except top border
            if ix == 0 or ix == nx-1 or iy == 0 or iy == ny-1 or iz == nz-1:
                density = 9999
            else:
                density = 2.5  # g/cm3
            row = {
                'ix': ix,
                'iy': iy,
                'iz': iz,
                'X': x_coord,
                'Y': y_coord,
                'Z': z_coord,
                'cu_pct': cu_pct,
                'UG': ug,
                'Category': category,
                'Below_Topo': flagged_below_topo,
                'Density_g_cm3': density
            }
            rows.append(row)


df_export = pd.DataFrame(rows)
df_export.to_csv("block_model_final_flagged.csv", index=False)

## for SGEMS quick validation
GSLIB.Dataframe2GSLIB("block_model_final_flagged.gslib", df_export)